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ABSTRACT 

Faraday Rotation Measure (RM) synthesis requires the recovery of the Faraday Dis- 
persion Function (FDF) from measurements restricted to hmited wavelength ranges, 
which is an ill-conditioned deconvolution problem. Here, we propose a novel decon- 
volution method based on an extension of the MUltiple Signal Classification (MU- 
SIC) algorithm. The complexity and speed of the method is determined by the eigen- 
decomposition of the covariance matrix of the observed polarizations. We show nu- 
merically that for high to moderate Signal to Noise (S/N) cases the RM-MUSIC 
method is able to recover the Faraday depth values of closely spaced pairs of thin 
RM components, even in situations where the peak response of the FDF is outside 
of the RM range between the two input RM components. This result is particularly 
important because the standard deconvolution approach based on RM-CLEAN fails 
systematically in such situations, due to its greedy mechanism used to extract the RM 
components. For low S/N situations, both the RM-MUSIC and RM-CLEAN methods 
provide similar results. 
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1 INTRODUCTION 

Faraday rotation is a physical phenomenon where the po- 
sition angle of linearly polarized radiation propagating 
through a magneto-ionic medium is rotated as a func- 
tion of frequency. Faraday Rotation Measure (RM) syn- 
thesis is an important method for analyzing multichan- 
nel polarized radio data, where multiple emitting regions 
are p resent along the single line of sight of the observa- 
tions (jBrentiens fc de Bruvnll2005l V In practice, the method 
requires the recovery of the Faraday Dispersion Function 
(FDF) from measurements restricted to limited wavelength 
ranges, which is an ill-conditioned deconvolution problem, 
raising important computational difficulties. At least four 
different approaches have been proposed so far to solve 
this problem. A firs t approach use s an adaptation of the 
CLEAN algorithm JHogboml Il974l) to the RM deconvo- 
lution (RM-CLEAN) (|Heald et al.lboogi '). The second ap- 
proach is wavelet-based, and assumes field sy mmetries in 
order to project the observed data onto A^ < (|Frick et al.l 
I2OI0I). The third approa ch relies on nonlinear model fitting 
( Farnsworth et al.lboill 'l , while the fourth approach is based 
on the novel compressed sensing (CS) paradigm (jLi et al.l 
l201ll : [Aridrecut et al.ll20ia '). Whether these methods are suc- 
cessful or not in detecting the structure of the FDF depends 
on the Signal to Noise ratio (S/N), the separation of the 
RM components, the relative angle of the RM components 
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at the observed wavelengths, and the wavelength range. For 
example, the current existing methods fail to recover the 
correct RM components when the separation among them 
is below the resolution determined from the half maximum 
of the main peak of the RM spread function (RMSF). Also, 
another major problem is related to situations where the in- 
terference of pairs of RM components conspire to place the 
peak response of the FDF outside of the RM range between 
the two input RM components. In this case, the standard 
RM-CLEAN method fails, due to its intrinsic greedy mecha- 
nism used to pick up the RM components. Here we discuss a 
novel approach, which addresses these two significant prob- 
lems from a different perspective. The proposed method is 
an extension of the MUltiple Signal Classification (MUSIC) 
algorithm, which is based on the eigen-deco mposition (ED ) 
of the covariance matrix of the observed data (|Chenevl200ll ). 
The complexity and the speed of the method is therefore 
determined by the size of the eigen-decomposition problem, 
which for several hundreds of data points is comparable to 
the RM-CLEAN approach. Our numerical results show that 
for high to moderate S/N cases, the RM-MUSIC method 
gives very good results, outperforming the standard ap- 
proach based on RM-CLEAN. For low S/N situations, both 
the RM-MUSIC and RM-CLEAN methods provide similar 
results. We should note that contrary to the other existing 
methods, the RM-MUSIC method recovers only the Fara- 
day depth values. Once the Faraday depth values are deter- 
mined, the real and imaginary parts of the RM components 
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can be easily computed using the linear least squares fitting 
approach. 



F{<P) = Y. f" 



(9) 



2 RM-SYNTHESIS 

The Faraday depth (in radm~^) is defined as: 



0.81 



n„B ■ dr, 



(1) 



source 



where Ue is the electron density (in cm~ ) , B is the mag- 
netic field (in /xG), and dr is the infinitesimal path length 
(in parsecs). We also define the complex polarization as: 



PiX") = Q{X') + iU{\^) = p/e^''^'^'\ 



(2) 



where p is the fractional polarization, /, Q, U are the ob- 
served Stokes parameters, and x(-^^) is the polarization an- 
gle observed at wavelength A. Also, we assume that the ob- 
served polarization -P(A^) originates from the emission at all 
possible values of <j), corresponding to the Fourier transform: 



P(A^) = 



F{(t>)e^"'^ d<f, 



(3) 



where F{(f)) is the complex FDF (the intrinsic polarized flux, 
as a function of the Faraday depth). Thus, in principle F{(j)) 
is the inverse Fourier transform of the observed quantity 



FW = 



P{\')e-^"^^' 



dy 



(4) 



However, this operation is ill-defined since we cannot observe 
P(A^) for A^ < 0, and also in general the observations are 
limited to an interval [\min , '^max] ■ 

In order to deal with the above limitations, the observed 
polarization is defined as: 

v2\ i!fr/\2\r}r\2\ 



p{x^) = w{y)p{x^), 



(5) 



where W is the observation window function, with VF(A^) > 
for A^ e [Xl,,„,Xl^ax], and W{X^) = otherwise. There- 
fore, for A'^ measurement channels with frequencies fn, wave- 
lengths An = c/fn, n — 1,2, ...,N (c is the speed of light), 
weights W{X^) = Wn, polarizations -P(A^) = P„, we obtain 
the following discrete expression for the FDF: 

F(0) = 



JV-l 
n = 



W„Pne 



-2i4,(\^-Xf.) 



(6) 



where 



A 



N-l 

E 



w„ 



(7) 



is a normalization constant, and the reference wavelength 
A^ is defined as: 



Considering additive noise in the measurement process, the 
sampled polarization-domain channel response is given by: 

K 

Pn = W„J2 /fce''^'=<''""^''-Ha„+i/3„, n = 1, 2, ..., N, (10) 

where A'' is the number of measured channels. Here, both a„ 
and Pn are Gaussian variables with mean zero and standard 
deviation a. Obviously, F{(j)) is a "dirty" reconstruction, and 
a deconvolution step is necessary to recover the RM com- 
ponents fkS{<l) — 4>k), given the observed values Pn and the 
weights Wn- 



3 RM-MUSIC 

The MUSIC method is generally used for estimating fre- 
quencies in signal processing problems, and the loca t ion o f 
pointlike scatterers in imaging problems fsee lChenevI (120011 ) 
for a review) . The standard MUSIC method is based on the 
ED of an Hermitian operator, which corresponds to the co- 
variance matrix of the signal or the array response matrix. 
By the finite-dimensional spectral theorem, such operators 
can be associated with an orthonormal basis of the under- 
lying space in which the operator is represented as a di- 
agonal matrix with real number entries. The main idea is 
to estimate the frequencies, or to localize multiple sources, 
by exploiting the eigen-structure of this Hermitian opera- 
tor. More exactly the space spanned by its eigenvectors can 
be partitioned into two orthogonal subspaces, namely the 
signal subspace and the noise subspace. By exploiting the 
orthogonality of the signal and noise subspaces, the MUSIC 
method significantly improves the resolution (i.e. locating 
closely spaced frequencies or scatterers), and as a conse- 
quence it is considered a super-resolution method. Let us 
now formulate the standard MUSIC approach to the Fara- 
day depth recovery problem. 

Since only one snapshot of measurement data P of 
length N is available, the data sequence is divided into 
M = N ~ L segments of length L, < L ^ N, and then the 
L X L covariance matrix is estimated as: 



Q^_Ly^p(-)(P 



{m)xH 



(11) 



where P'"' = [Pm,Pm+i, -, Pm+L^ , m = 0,1,...,M- 1, 
and the superscript H denotes the conjugate transpose 
operation. Suppose that fJ,o ^ fJ,i ^ ■■■ ^ ^J■L-l and 
g'° , g , ..., g'^~^' are the eivenvalues, and respectively the 
eigenvectors of G, such that: 

and 



(12) 



A? = A-'^vy„A^. 



(8) 



In our approach we assume that the model of P(<jf>) con- 
tains K <^ N (unknown) components fkS{(ji~ (pk), where fk 
are complex and (pk are real quantities fc = 0, ..., A" — 1: 



£=0 



(13) 



Since the measured signal contains only K components, the 
last L — K eigenvalues of G, fj.K ^ I^-k+i ^ ... ^ I^-l-i, 
should be small and below the noise level, and we say that 
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the corresponding eigenvectors g'^^^g'^^^', ..., g'^^^' span 
the noise subspace of G, while the first K eigenvectors 
g' , g' ,...,g' ~ •* span the signal subspace, and the cor- 
responding eigenvalues and eigenvalues /io ^ /ii ^ ... ^ 
/iif-i should be above the noise level. Thus, the eigen- 
decomposition of G can be written as: 



G=^M.g<^'(g'^')^ + X,M.g^^^(g^^^)^ 



rW(„W\'^ 



(14) 



signal subspace 



noise subspace 



We can now form the projection operator onto the noise 
subspace defined as: 



G 



{noise) 



E 



s''Hk''Y, 



and consider the signal subspace sampling vector 

h(0) = [e^'-^^o^e'^'^^?, 



2i4>\i 



where the wavelengths are interpolated as following: 

-\ 2 



A 



UN - TTTii'^N - Vl) 



0,1, 



1, 



(15) 



(16) 



(17) 



such that they cover the initial observation interval 
[\min,^max\ We should uotc that since L ^ N, the Fara- 
day depth is also linearly interpolated as </) = ij>L/N , and 
therefore we obtain: 



h(0) = [e "'^0,6 N^ 1 ^ 



(18) 



If the vector h(0) represents the signal, i.e. it is a linear 
combination of the signal subspace eigenvectors, then its 
projection onto the noise subspace must be close to zero: 



i {noise) i 



|G^"°"'=^h(0) 



0. 



(19) 

Thus, the Faraday depth values (j)k, k = Q,1, ..., K—l, should 
correspond to the maxima of the following MUSIC pseudo- 
spectrum: 

1 1 



S{4>)^ 



J^ (noise) 



g 



h(0) 



E 



''-'\\isW)Hhm'' 



(20) 



Once the Faraday depth values (j)k are determined from 
the MUSIC pseudo-spectrum the real and imaginary parts 
of the RM components fk can be easily computed using 
the linear least squares fitting approach, i.e. by solving the 
following minimization problem: 

-I 2 



{/fc}fc=o =argnim 



fk 



iV-l 

E 



K 

E 



/fee 



'2i't>k^„ 



= arg min P 
f " 



*f 



(21) 



This least-squares fit problem has fewer free parameters, and 
it is far better constrained than the fit that would have been 
done without MUSIC. In fact the problem has an unique 
minimum norm solution, obtained by solving the linear sys- 
tem of equations: ^f = P, where f — [/o, /i, ..., Jk-i]'^ are 
the unknown components, P = [Po, Pi, ..., Pjv_i]'^ are the 
observed polarizations, and \E' is the Fourier matrix with 
the elements '$n,k = e^"^'='^". Thus, the least squares solu- 
tion is given by f = *tp^ where *+ = (*^*)-i*^ is the 
Moore- Penrose pseudo-inverse of SE*. 



Algorithm 1 Pseudo code of the RM-MUSIC method. 

P; polarization data 

A^; number of channels 

N/2 ^ L ^ N; length of data segments 

M = N — L; number of data segments 

G ^ 10- EmJo P(™'(P<'"')^; covariance matrix 






ED{G); eigen-decomposition 

1 2 



VN-j^-ri('^N-'^l) 



h(0) 



„2i^0A^ g2i^0A^ 



; interpolation wavelengths 

p 'Tv"* i-i]-^; sampling vector 



K -ir- arg minx H{K); number of components 



S{<t>) 



MUSIC pseudo-spectrum 



{(t>k}kJo ^ {argmax^S(cj!>)}f^(,^; Faraday depths 

P — '^f ; complex amplitudes 

\K-1. 



{/fejf^o^ ^argmiuf I 
return {(pk,fk}kJo 



4 NUMERICAL IMPLEMENTATION 

In order to optimally use the MUSIC method we also need 
to determine the number of components K in the polar- 
ization signal. In principle, K can be determined from the 
eigenvalues of the covariance matrix. The first largest eigen- 
values, which are much bigger than the noise level /ii > a^, 
i? = 0, 1, ...K — 1, should indicate the value of K. However, 
in a practical implementation, when the covariance matrix 
is estimated from a small number of observations, it is chal- 
lenging to clearly separate the signal eigenvalues from the 
noise eigenvalues. In these ambiguous cases one can use in- 
formation theoretic criteria for selection, like the Minimum 
Descriptive Length (MDL) criteria, where the value of K 
corresponds to the minimum of the following quantity: 



H{K) = -M{M - K) log 



+ -K{2L- K)logM. 



YfM 1/(M-K) 



l-K L-i 



M-K L-im=K 



{-I'm 



(22) 



Also, since the smallest eigenvalues are not equal among 
them, one can use the following normalized versions of the 
MUSIC pseudo-spectrum: 



5(0) 



EtKM7Mi(gW)^h(<^)ir 



(23) 



that accounts for the variation of the eigenvalues, and it is 
less sensitive to the K estimation errors. Another infiuencing 
parameter is the length of the data segments L. In order to 
preserve as much information as possible in each data seg- 
ment we should have L ^ A'^/2. Our numerical experiments 
have shown that L = 2N/3 has a good detection sensitiv- 
ity, and this is the value used in the simulations presented 
here. We should note also that the speed of the RM-MUSIC 
depends on the parameter L, which gives the size L x L of 
the eigen-decomposition problem. For L in the order of hun- 
dreds, the speed of RM-MUSIC is comparable to the speed 
of RM-CLEAN, or even faster. The pseudo-code of the RM- 
MUSIC method is given in Algorithm 1. 
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Figure 1. RM-MUSIC successful result for two components 
with the complex amplitudes /i = — lOOi, /i = 50i, and Fara- 
day depths (po = Oradm"^, 0i = 122radm~^, separated by 
4i\ — 4'2 — <5</> (I/I component amplitudes, |i^| "dirty" FDF, |c| 
RM-CLEAN components, S RM-MUSIC spectrum). 



Figure 2. RM-MUSIC successful result for two components 
with the complex amplitudes /i = — lOOJ, /i = 50J, and Fara- 
day depths (/iQ = Oradm"^, ij>\ = 61radm~^, separated by 
</>! — 02 — 0.55(/) (I/I component amplitudes, |i<'| "dirty" FDF, 
c| RM-CLEAN components, S RM-MUSIC spectrum). 



5 NUMERICAL RESULTS 

To illustrate numerically the RM-MUSIC method we con- 
sider the following experiment layout: frequency range: 
I'min = 1100 MHz, Umax ~ 1400 MHz; wave length range 



AL„ = 0.046 m^, 



^max = 0.074 m^; number of observa- 



tion channels A'^ = 150; width of an observing channel 
Av = 2 MHz; resolution in Faraday depth space 5(j) = 
^Vs/^Xmax — A^i„) ~ 122radm~^; equal weights W„ — 1; 
noise level a = \/F = 12.247. 

The RM-MUSIC method works very well when the com- 
ponents are separated above the resolution limit 5(j). For ex- 
ample, in Figure 1 we have two components with the com- 
plex amplitudes /o = — lOOi and /i = 50i at A^ — A^, 
and Faraday depths (j>o = Oradm"^, 01 — 122radm~^, 
separated by cpi — 4>o = 5(j). For comparison we also give 
the "dirty" FDF and the RM-CLEAN components. In or- 
der to make the results comparable, we scaled S{(f)) as 
S(</.) ^ S(0) |F|_^ |5|-L, where \F\^^^ and |S|_, are 
the maximum amplitudes of the original spectra. This way, 
both S{'j)) and 1-^(0)1 are in the same range. One can see 
that the RM-MUSIC spectrum recovers almost exactly the 
Faraday depths of the two components (^o = 1.05radm~^, 
4>i = 122.81 radm"^ while both the "dirty" FDF and RM- 
CLEAN fail to recover the Faraday depth of the second com- 



ponent: 



2.01 rad m" 



142.95 radm .Moreover, 



once the depths are recovered, the real and imaginary parts 
of the RM components can be easily computed using a linear 
least squares fitting approach, resulting in an almost exact 
recovery: fa = -1.95 - 99.74i, /i = 0.05 + 50.02i. 

In Figure 2 we show that the MUSIC super-resolution 
method can give good results even for separations be- 
low the the resolution in Faraday depth space S(f>. Here, 
we have the same two components with Faraday depths 
(f)o = Oradm"^, (f>i = 61radm~^, separated by (j>i — cj>o = 
0.550. One can see that both components are correctly 
separated by the RM-MUSIC spectrum 0o = Oradm"'^, 
/o = 0.84-100.93i, 4>i =62.4radm-^ /i = -0.76-h51.49i. 
and once again the "dirty" FDF and RM-CLEAN fail to 
recover the Faraday depths, showing two main components 
at: 00 = — 6.4radm~^, and respectively 0i — 118.8 radm~^. 
This is a typical example where the interference of pairs of 
RM components conspire to place the peak response of the 
FDF outside of the RM range between the two input RM 
components. In Figure 3 we consider the same two com- 
ponents separated by 0i — 0o — 0.25(50, 0o ~ Oradm"^, 
01 — 30.5 rad m~^. In this case the RM-MUSIC spectrum 
cannot separate the components anymore, however it cor- 
rectly shows only one majcimum, situated between the two 
input RM components, while RM-CLEAN suggests a single 
RM component outside the range of the input RMs. 

From the above examples we see that RM-MUSIC per- 
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Figure 3. RM-MUSIC result for two components with tiie com- 
plex amplitudes /i = — lOOi, /i = 50i, and Faraday depths (f>o = 
Oradm"^, 0i = 30.5radm~^, separated by <?ii — </>2 — 0.255<j) (|/| 
component amplitudes, \F\ "dirty" FDF, |c| RM-CLEAN compo- 
nents, S RM-MUSIC spectrum). 
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Figure 4. RM-MUSIC result for low signal to noise is compara- 
ble to RM-CLEAN. Two components with /i = -lOOi, /i = 50i, 
Faraday depths 0o = Oradm"^, (f>i = 61radm~^, separated by 
4>i—4>2— 0.b&4> (I/I component amplitudes, \F\ "dirty" FDF, |c| 
RM-CLEAN components, S RM-MUSIC spectrum). 



forms better than RM-CLEAN for high to moderate S/N 
cases. However, in low S/N situations the performance of 
RM-MUSIC deteriorates. It is quite difficult to quantify the 
S/N ratio exactly and to find a limiting S/N value, and more 
future evaluations of the method are required for various 
wavelength ranges and noise levels per channel. For the con- 
sidered experiment layout our numerical simulations have 
shown that if the amplitude of the strongest RM compo- 
nent is becoming smaller than |/oj ~ 50, and for a noise 
per channel a = y/N = 12.247, RM-MUSIC begins to 
behave more like RM-CLEAN, and cannot separate cor- 
rectly the input components. In Figure 4 we give such an 
example of two components with the complex amplitudes 
/o = — 30i and /i — 15i, Faraday depths ^o = radm~^ and 
4>\ = 61 rad m^^, separated by 4>i — 4'o — 0.5S4>. One can see 
that in this case, both RM-MUSIC and RM-CLEAN give 
similar results. We should also mention that the experiment 
layout used in si mulations is very sim ilar to the ASKAP's 
POSSUM survey (JGaensler et aLlboid ). where the expected 
final noise is 10 /ijy/beam. This value suggests that for 
A'^ = 150 channels and a peak signal to noise threshold of 
~ 4, RM-MUSIC should work as indicated for sources with 
polarized flux higher than 0.5 mjy/beam, which is rather 
bright but quite reasonable, since about ~ 57% of the NVSS 
RM catalog sources fall into this category. 



6 CONCLUSION 

We have discussed a novel deconvolution method for RM 
synthesis applications, based on an extension of the MUSIC 
super-resolution algorithm. Numerical results show that for 
high to moderate S/N cases RM-MUSIC outperforms the 
standard RM-CLEAN approach, being able to recover Fara- 
day depth values of closely spaced pairs of thin RM compo- 
nents, even in situations where the peak response of FDF is 
outside the range between the two input components. For 
low S/N cases, the RM-MUSIC and RM-CLEAN methods 
provide similar results. RM-MUSIC performs well for pairs 
of thin RM components, however it is not suitable when 
the input Faraday spectrum contains polarized emission at 
a continuous range of RM. 
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